Strongly correlated metal interfaces in the Gutzwiller approximation 
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We study the effect of spatial inhomogeneity on the physics of a strongly correlated electron 
system exhibiting a metallic phase and a Mott insulating phase, represented by the simple Hubbard 
model. In three dimensions, we consider various geometries, including vacuum-metal-vacuum, a 
junction between a weakly and a strongly correlated metal, and finally the double junctions metal- 
Mott insulator-metal and metal-strongly correlated metal- metal. We applied to these problems the 
self-consistent Gutzwiller technique recently developed in our group, whose approximate nature is 
compensated by an extreme fiexibility, ability to treat very large systems, and physical transparency. 
The main general result is a clear characterization of the position dependent metallic quasiparticle 
spectral weight. Its behavior at interfaces reveals the ubiquitous presence of exponential decays and 
crossovers, with decay lengths of clear physical significance. The decay length of metallic strength in 
a weakly-strongly correlated metal interface is due to poor screening in the strongly correlated side. 
The decay length of metallic strength from a metal into a Mott insulator (or into vacuum) is due to 
tunneling. In both cases, the decay length is a bulk property, and diverges with a critical exponent 
(-^ 1/2 in the present approximation, mean field in character) as the (continuous, paramagnetic) 
Mott transition is approached. 



I. INTRODUCTION 

Metallic electron wavefunction delocalization in a lat- 
tice of atoms or molecules is caused by the lowering of 
electron kinetic energy and by the simultaneous improve- 
ment of electron-ion Coulomb attraction. By abandon- 
ing the ion cores and turning delocalized, an electron 
can in fact feel the potential of more than one nucleus. 
However, coherent electron motion is opposed by the mu- 
tual electron-electron Coulomb repulsion, which is higher 
when electrons move due to their higher chance of col- 
liding when visiting the same site. When the first two 
terms prevail, the system is a conventional band insu- 
lator or metal, depending whether the Fermi level falls 
in a band gap or across one or more bands. When the 
electron-electron repulsion prevails instead the electrons 
localize on their atomic or molecular sites leading to a 
so-called Mott insulator—. Despite that conceptual sim- 
plicity, properties of Mott insulators and especially of 
strongly correlated metals in the proximity of a Mott 
metal-insulator transition as a function of increasing cor- 
relations remain quite difficult to capture both theo- 
retically and experimentally. Theoretically, the reason 
is that the Mott transition is a collective phenomenon, 
which escapes single-particle or mean field theories such 
as Hartree-Fock or density-functional-theory within the 
local-density approximation (LDA). Experimentally, ad- 
ditional complications such as magnetism, lattice distor- 
tions, etc., often conspire to mask the real nature of the 
Mott localization phenomenon. 

Important insights into this problem have been gained 
in the last two decades especially thanks to dynamical 
mean field theory (DMFT),^ DMFT predicts that, as 
the electron-electron repulsion - usually parametrized by 
a short-range Hubbard repulsion U - increases, the or- 
dinary band metal evolves first to a strongly correlated 



metal well before the Mott transition. In the strongly 
correlated metal the electron spectral function under- 
goes a profound change, exhibiting well formed Mott- 
Hubbard side-bands coexisting with delocalized quasi- 
particles, the latter narrowly centered in energy near the 
Fermi level. Only successively upon increasing repulsion 
do the quasiparticles disappear as the Mott transition 
takes place at C/ = Ud-it- This intriguing prediction - 
simultaneous metallic and insulating features, exhibited 
on well separated energy scales - has stimulated a con- 
siderable experimental effort to reveal coexisting quasi- 
particles and Mott-Hubbard bands in strongly correlated 
metals^'^'^'i'^'^'iSiiii'i^, especially in the paradigmatic 
system V2O3. This is the compound where a Mott tran- 
sition has been first discoveredp^ and theoretically stud- 
Ig ji4,i5 _ ambient temperature and pressure V2O3 is a 
correlated metal. It undergoes a first-order Mott transi- 
tion at ^ T/v ~ 155 K to an antiferromagnetic insulator 
accompanied by a monoclinic distortion of the high tem- 
perature corundum structure.— The paramagnetic high- 
temperature metal can moreover be turned into a param- 
agnetic Mott insulator upon substituting V with bulkier 
Cr, (Vi_:,Cr3,)203. For 0.005 < a; < 0.017 a first-order 
line separates the high temperature metal from the para- 
magnetic Mott insulator, which terminates with a critical 
point at T ~ 400 K and x ~ 0.005. 

Near the metal-insulator transition of (Vi_a;Cr2:)203, 
the strongly correlated metal must of course possess well 
defined quasiparticles at the Fermi energy. Surprisingly, 
early photoemission experiment a^^i^^i^^i^° failed to re- 
veal the sharp quasiparticle peak predicted by DMFT 
at Ef- The electronic spectrum appeared instead dom- 
inated by the lower Mott-Hubbard band with barely a 
hint of metallic weight at the Fermi energy. It was rec- 
ognized only later that photoemission in strongly corre- 
lated metals is highly surface-sensitive j^i^il i 1 1 1 1 2 , 2 1 gy 
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increasing the photon frequency, which corresponds to 
more energetic excited photo-electrons, i.e. longer es- 
cape lengths, a prominent quasiparticle peak coexist- 
ing with incoherent Mott-Hubbard bands was eventu- 
ally observed in ¥203^ '^°'^^ . Quasiparticle suppression 
in surface-sensitive probes was attributed^^ to surface- 
modified Hamiltonian parameters, the reduced atomic 
coordination pushing the surface closer to the Mott tran- 
sition than the underlying bulk. This conclusion, al- 
though not unreasonable, raises however a more funda- 
mental question. A metal does not possess any intrinsic 
long-distance electronic length-scale other than the Fermi 
wavelength. Thus an imperfection like a surface can only 
induce a power-law decaying disturbance such as that 
associated with Priedel's oscillations. Since one does not 
expect Luttinger's theorem to break down, these oscilla- 
tions should be controlled by the same Fermi wavelength 
as in the absence of interaction, irrespectively of the 
proximity of the Mott transition. On the other hand, a 
strongly correlated metal does possess an intrinsic energy 
scale, the parametric distance of the Hamiltonian from 
the Mott transition, and that could be associated with a 
length scale. For example, the arising of a critical length 
scale in association with a free energy scale is well known 
in second order phase transitions. The surface as a per- 
turbation may alter the quasiparticle properties within 
a depth corresponding to that characteristic length. We 
expect this length to be a bulk property, the longer the 
closer the Mott transition, unlike the Fermi wavelength 
that remains constant. In this respect, it is not a priori 
clear whether the recovery of bulk-quasiparticle spectral 
properties with increasing depth should be power-law, 
compatible with the common view of a metal as an in- 
herently critical state of matter, or exponential, as one 
would expect by regarding the Mott transition as any 
other critical phenomena where power laws emerge only 
at criticality. 

Besides the interface with vacuum, which is relevant 
to spectroscopy, other types of interface involving cor- 
related materials are attracting increasing interest. In 
2004, Ohtomo and Hwang^^ discovered that the inter- 
face between two insulating oxides, LaAlOa and SrTiOa, 
is a high-mobility two-dimensional conductor that even 
shows superconductivity^!. This discovery stimulated ex- 
perimental and theoretical studies on oxides heterostruc- 
tures^^. On the theory side, some activity has been 
focused either on the characterization of the electronic 
structure of these interfaces by ab-initio LDA calcula- 
tions, see e.g. Ref. [H, as well as on DMFT analy- 
ses of simple model921i2£i2£i^i^i2^iHiMi^ and on com- 
bined LDA-DMFT calculations^ aimed at understand- 
ing interface correlation effects poorly described within 
straight LDA. The DMFT approaches adopted in the lit- 
erature to describe this kind of situations were ad-hoc ex- 
tensions of the single-site DMFT— to inhomogeneous sys- 
tems.— In the specific example of a layered structure, 
the electron self-energy was assumed to depend, besides 
the frequency, also upon the layer index. In this scheme 



the self-energy is calculated by solving an auxiliary impu- 
rity model for each layer in which the conducting bath de- 
pends self-consistently on the fully-interacting impurity 
Green's functions not only of that given layer but also 
of the nearby ones. This additional complication with 
respect to conventional DMFT weighs on the numerical 
calculation, which is thus limited to few tens of layers. 
Although this is adequate for the interface between two 
insulators, such as that studied by Ohtomo and Hwang^, 
it is generally insufficient in other cases, such as the sur- 
face effects in the interior of a correlated metal,— or any 
other interface involving at least one metal. 

Recently, we proposed an alternative theoretical ap- 
proach to interface problems,— based on the extension 
of the Gutzwiller wavefunction and approximation^^iiS 
to inhomogeneous situations. The method, although a 
further approximation beyond DMFT, hence in principle 
less accurate, is much more agile, and can treat without 
effort hundreds of layers. Thus it can be used as a com- 
plementary tool to extrapolate DMFT results to large 
sizes, otherwise unaccessible by straight DMFT. 

In this work, we shall extend the analysis of Ref. [H for 
the vacuum / correlated- metal interface to other model in- 
terfaces that might be relevant for experiments: the junc- 
tion between two different correlated metals and the tun- 
neling between two metallic leads through a strongly cor- 
related, possibly Mott insulating, region. Although both 
cases were in fact previously studied by DMFT— i^i^, 
the results were interpreted in contrasting ways. While 
Helmes et ai.— concluded that the Mott insulator is im- 
penetrable to the electrons coming from the metallic 
leads, Zenia et al^ drew the opposite conclusion that 
a conducting channel always open up inside the insula- 
tor at sufficiently low temperature. The present study, 
which is certainly less accurate than DMFT but can deal 
with much larger sizes, will also serve to clarify this is- 
sue. In particular, the large sizes allow us to address the 
asymptotic behavior and to identify the magnitude and 
interface role of the critical length associated with the 
bulk Mott transition. 

The paper is organized as follows. In section lUl we 
introduce the model Hamiltonian, which is a Hubbard 
model with layer dependent parameters, and a Gutzwiller 
variational scheme adapted for such an inhomogeneous 
situation. We then study in sectionlUllthree different slab 
geometries: (a) strong correlated metal-vacuum inter- 
face; (b) junction between two different correlated met- 
als; (c) a Mott insulator or a strongly correlated metal 
sandwiched between two weakly correlated metals. In 
the first two cases we find that the perturbation induced 
by the surface inside the bulk of the correlated metal de- 
cay exponentially at long distances. The length scale ^ 
that controls this decay is a bulk property that depends 
in our simplified model only on Ucrit — U and diverges on 
approaching the Mott transition like ^ ^ (?7crit — U) , 
with a mean-field like exponent i' ~ 0.5. The last case 
(c) is more interesting. Either when the central region, 
of width d, is a strongly correlated metal, C/contcr < C^rit 
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or when it is a Mott insulator, L'contor > UcHt, the effects 
of the two metal leads are found to decay exponentially 
over a length ^. Just like in cases (a) and (b) above, ^ is 
only controlled by the distance from Mott criticality, i.e. 
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which therefore appears naturally as a correlation length 
that is finite on both sides of the transition. However, 
while the quasiparticle weight saturates to a finite con- 
stant determined by J7center < K:rit and independent of 
d when the central region is a strongly correlated metal, 
in the opposite case of a Mott insulator the quasiparti- 
cle weight saturates to a finite value exponentially small 
in d. Interestingly, right at criticality, C/contcr = f^crit, 
the saturation value decays power law in d. Finally, sec- 
tion Uv] is devoted to concluding remarks. For a better 
understanding of our numerical data, a simple analytical 
model for the spatial dependence of quasiparticle weight 
is set up in appendix O while in appendix [B] we discuss 
the effects of electron-electron interaction on the physics 
of Friedel's oscillations near surfaces and junctions within 
the Gutzwiller approximation. 



II. MODEL AND METHOD 

In order to address the generic interface features of 
a a strongly correlated metal, we consider the simplest 
Hamiltonian exhibiting a Mott transition, namely the 
Hubbard model 

<RR'>cr 

+ ^ CRflR + Ur riRT'^Ri, (1) 
R 

where < RR' > denotes nearest neighbor sites, c^^ 
and Cj^^ creates and annihilates, respectively, an elec- 
tron at site R with spin a, and finally n-Ra ~ c^o-'^Ro- 
and riR — nR| -f- riRj . In our inhomogeneous system, 
all Hamiltonian parameters are allowed to be site depen- 
dent. For interfaces, we shall assume an iV-layer slab 
geometry where all parameters are constant within each 
layer, identified by a layer coordinate z = 1, . . . ,N but 
generally different from layer to layer. For instance, the 
hopping between nearest neighbor sites R and R' within 
layer z depends only on z, i.e. tRR' = t{z), while if R 
and R' belong to nearby layers, e.g. z and z ± 1, then 

tRR' =t{z,Z±l)^t{z±l,z). 

We study the Hubbard Hamiltonian (U) in the non- 
magnetic (also called paramagnetic) sector by means of 
a Gutzwiller type variational wavefunction 



I*) =n^ai^o), 



(2) 



R 



where |^'o) is a paramagnetic Slater determinant. Be- 
cause of our choice of layer-dependent parameters, the 



operator Vr has the general expression 
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Vr^Y. A„(z)|n,R)(7i,R|, 



(3) 



n=0 



where |7i, R)(n,R| is the projector at site R = {x,y,z), 
{x and y are intralayer coordinates), onto configu- 
rations with n electrons (note that |1,R)(1,R| = 
c^^|0, R)(0, R|cj^^), and A„(z) are layer-dependent 
variational parameters. We calculate quantum averages 
on \^) using the so-called Gutzwiller approximatio n^^i*^" , 
(for details see e.g. Ref4i whose notations we use here- 
after) and require that 



(*o|7'^|*o) 
(*o|:PRnR.|*o) 



(*o|"R<t|*o) 



n{z) 



(4) 
(5) 



Explicitly, these two conditions imply that 
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1 



n(z) 



n[z) 



+n{z)[l-^\x,{zf + '^\,{z)\{(S) 



(7) 



niz)[l-mX,izf + 2^-^Uzf. 



We note that n{z) is fixed once the uncorrelated vari- 
ational wavefunction |5'o) is given. In reality we find 
more convenient to treat n{z) as an additional variational 
parameter, and constrain |^'o) to span all paramagnetic 
Slater determinants that have a fixed local charge den- 
sity n{z). The average value of (U) within the Gutzwiller 
approximation is accordingly given by^>^ 

^ = -^^''T.UT.^\2{zf + eRn{z) (8) 

R. 

- iRR'i?(^)i?(^')(*o|cLcR-. + i?.c.|«'o), 

<RR'><T 



""^^^^ Ao(z)Ai(z) + ^Ai(z)A2(z), (9) 



where 

R{z) - ( 1 



plays the role of a wavefunction renormalization fac- 
tor, whose square can be regarded as the actual layer- 
dependent quasiparticle weight, Z{z) = R^{z). Because 
of Eqs. ([6]), ^ and ([9]), one can express 

A„(z) = A„ [R{z),n{z)], 

as functional of the two variational functions R{z) and 
n{z). Furthermore, the single-particle wavefunctions that 
define the Slater determinant |5'o) can be chosen, for a 
slab geometry, to have the general expression 



'?!'ek||(R) 



iki I -R, 



'/'ekM(^;), 
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where A is the number of sites per layer and k| | the mo- 
mentum in the x-y plane. The minimum of E, Eq. l(8|. 
can then be obtained by searching for saddle points with 
respect to the variational parameters R{z), n{z) and 
0£k|| {z), the latter subject to the constraint 

2 occupied 

^ \4'ek^^{z)\ =n{z), 

the sum running over all occupied states in the Slater 
determinant. 



Considerable simplifications arise if we further assume 
a bipartite lattice with a Hamiltonian ^ invariant under 
the particle-hole transformation 

where (—1)''^ is -1-1 on one sublattice and —1 on the other. 
This symmetry requires eR = in ((!]) and implies n{z) = 
1 hence Ao(z) = ^2(2;) and Xi{z)'^ = 2 — Ao(z)-^. In this 
case the saddle point is simply obtained by solving the 
coupled equations 



e0ek||(2) = i?(z)^ek||(2)0£k||(2:) - Riz) ^ t{z, z + p) R{z + pa) (f)ei^^^{z + pa), 

P=±l 



A 777 — To occupied 

Rjz) = '^]^ E 



U{z)A 



eki 



- 2i?(z)ek||(2)</'ek||(z)^ -I- 0£k|| (-z) E t{z,z+p)R{z+pa)(j)f_\^^^{z+pa) 

p—±a 



(10) 

,(11) 



where ek||(2) = —2t{z) {cos k^a + cos kya) . The first 
equation has the form of a Schroedinger equation which 
the single-particle wavefunctions ^ek, | (z) must satisfy, 
the quasiparticle hopping now depending parametri- 
cally on R{z). The second equation has been in- 
tentionally cast in the form of a map Rj+i{z) = 
F [Rj{z), Rj{z + a), Rj{z — a)] whose fixed point we have 
verified to coincide with the actual solution of (fTTI) in the 
parameter region of interest. 

In spite of the various assumptions above, solving this 
saddle point problem remains in principle formidable. 
Fortunately, Eqs. (flOl) and ifTTj) can in fact be solved rel- 
atively easily, by the following iterative procedure. First 
solve the Schroedinger equation at fixed Rj{z)\ next find 
the new i?j+i(z) using the old Rj{z) and the newly de- 
termined wavefunctions ^ekul-z)- With the new i?j+i(z), 
repeat the above steps and iterate until some desired level 
of convergence is reached. Because of the large number 
of variational parameters, this iterative scheme is much 
more efficient than - but fully equivalent to - a direct 
minimization of i?, Eq. |[8]). Away from particle-hole sym- 
metry, the saddle point equations get more involved but 
the solution can be obtained along the same lines. 

Before concluding, we recall for future use the 
Gutzwiller approximation results for the Mott transi- 
tion at particle-hole symmetry in the homogeneous case, 
ER = 0, iRR' = t and C/r = C/, i.e. when the varia- 
tional parameters A„(z) are ^-independent. In this case, 
the solution of Eqs. I|10p and (fTTj) is trivial. The critical 
values U = J7crit at the Mott transition are Ucmt = 32t/7r 
(for a linear chain), Ucnt — 128t/7r^ (for a square lattice), 
C/crit = 16t (for a cubic lattice). The quasiparticle weight 
Z in terms of the electron-electron interaction U has the 



simple expression 

Z = i?2 = 1 - , (12) 
linearly vanishing at the Mott transition.— 

III. INTERFACES IN THE 3D HUBBARD 
MODEL: RESULTS 

We use the technique just exposed to study 3D sim- 
ple cubic Hubbard model interfaces in a slab geometry 
with in-plane (xy) translational symmetry and layer(z)- 
dependent Hamiltonian parameters. We assume for sim- 
plicity particle-hole symmetry and site-independent hop- 
pings tRR' = t throughout, so that the only source of 
inhomogeneity is a layer-dependent U{z). Therefore the 
minimization procedure amounts to solve the coupled 
equations lfTO| and (fTTj) with constant hoppings. Tech- 
nically, we diagonalized the in-plane /c-dependent Hamil- 
tonian lfTO|) at every point of a Monkhorst-Pack fc-grid^. 
The two-dimensional grid used was 32 x 32, chosen so as 
to yield well converged values not just for the quasiparti- 
cle weight (for which a 4 x 4 grid was sufficient) but also 
for the hopping matrix element for the geometries and 
interaction parameters considered. At every iteration j, 
we choose for the convergence indicator 

a threshold of 10~^. This corresponds to a relative energy 
convergence of less than 10"''. The calculations of the 
spatial dependence of the hopping matrix elements (see 
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appendix[B]) were instead performed with a denser fc-grid 
of 64 X 64 fc-points. 

We consider the three different geometries displayed in 
Fig. [11 

(a) Correlated metal-vacuum interface: a correlated 
metal (C/buik < C^rit, where J/crit = 16i is the 
critical value of U at the Mott transition in the 
cubic lattice) with a stronger correlated surface 

(f^surface ^ 

(b) Weakly correlated metal-strongly correlated metal 
interface: a junction between a moderately corre- 
lated metal (C/ioft < t^crit) and a strongly correlated 
metal (C/right < C/crit)- 

(c) Metal-Mott insulator-metal double junction: a 
Mott insulator t/ccntcr ~ C4rit or a strongly cor- 
related metal J/ccntcr ^ C4rit sandwiched between 
two moderately correlated metallic leads U\eit = 

bright < C^rit- 

The dashed lines in the panels of Fig.[T]show the quasi- 
particle weight Z(z) calculated for a N = 200 layer slab 
in the three geometries with the Hamiltonian parameters: 

(a) C/buik - U{z > 1) = 15.9712t and C/,urfacc = U{z = 
1) — 20t. The bulk is a strongly correlated metal 
very close to the Mott transition, the right surface 
has the same U as the bulk while the left surface a 
higher value well inside the Mott insulating range. 

(b) Uictt = U{z < 100) = 15.9198t and C/right = U{z > 
100) = 15.9712i; The left metal is much less corre- 
lated than the right metal. 

(c) C/right = Uiz < 80) = C/ieft = Uiz > 120) = 
15.9198t and C/co„tcr = C/(80 < z < 120) = 
16.0288t. Left and right leads are moderately corre- 
lated metals, the central region is Mott insulating. 

We now discuss each case separately. 



A. Geometry (a): Correlated metal- vacuum 
interface 

This is the simple surface case, U{z > 1) = C/buik < 
C^crit and U{z = 1) = C/sm-faco > C/crit, prcviously studied 
in Ref. [H. Looking at Figs. [2] and [H with values of 
Usuri = 20i, and C/buik = 9-6t and C/buik = 15.97118i, 
respectively, we observe that: 

i) The value of Z{z) at the center of the slab, close 
to the bulk value, decreases monotonically to zero 
while C/buik approaches C/crit. Due to the finite 
slab thickness N, the actual value of U at which 
Z{z) vanishes everywhere is slightly smaller than 
the bulk value C/crit = 16i for an infinite system, 
but tends to it as N increases. In this Hmit, the 
dependence of ^buik = Z{z = N/2) upon C/buik is 
described by Eq. lfn]) . 
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Figure 1: (Color online) The three different inhomogeneities 
studied in this paper: (a) free surface geometry, (b) junc- 
tion between metals with different strength of correlation, 
(c) Mott (or strongly correlated metallic) slab sandwiched 
between metallic leads (sandwich geometry). The values 
for U in all the three cases shown cir6^ (cl) U surface — 20i , 

C/buik = 15.9712t; (b) !7ieft = 15.9198t, C/nght = 15.9712t; (c) 
(7ieft = C/right = 15.9198*, C/contor = 16.0288t (which is the 
case of a Mott central slab). In panel (c) the region with 
electron-electron interaction U = C/contcr is indicated by the 
green-shaded area. 
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Figure 2: Spatial dependence of Z{z) for C/surf ~ 20f at z = 
and C/buik ~ 14.6642*, for any z > 0. The lower panel is the 
same as the upper one zoomed close to the surface. 



a) Z{z) decreases dramatically while approaching the 
surfaces, both the extra-correlated left surface z = 
1, and the regular bulk-like one ai z = N . In fact, 
within the Gutzwiller approximation, the effective 
interaction strength at a given site is the value of 
U relative to the average hopping energy at that 
site. The reduced surface coordination lowers the 
overall hopping energy of a surface site, and hence 
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Figure 3: Same as Fig. [2] for (7suif = 20t and (7buik ~ 
15.9712t. 



effectively strengthens the surface interaction. The 
same effect would be obtained by decreasing the 
hopping at the surface. We note however that, so 
long as Z remains finite in the interior of the slab, 
Z remains finite, even if very small, also at the 
surface: there cannot be truly insulating surfaces 
coexisting with a metallic bulk. The reason is that, 
if we assume initially such an insulating surface, 
then simple tunneling from the underlying bulk will 
bring the metallic quasiparticle weight to a nonzero 
value, however small. 

Ill) The steep decay of Z{z) at the surfaces at 2 = 1 
and z — N gets more and more gradual as C/buik ^ 

As found in Ref. HI, the behavior oi R{z) ~ \J Z(z) 
can be well described by an exponential 

= i?bulk + (i?surf - i?bulk) e-("-l)/«, (14) 

where i?buik = ^(2 = ^/2) and i?surf < ^buik- In Ap- 
pendix A we actually derive a more involved analytical 
expression for R(z) that fits well the numerical data, see 
Eq. (jA6p . The surface value, i?surf , and the surface metal- 
lic quasiparticle weight ^surf = -Rsuifi ^re much smaller 
than the bulk ones but, as previously mentioned, they 
can vanish only when i?buik becomes strictly zero, for 
t/buik > C^crit- For any C/buik < C^crit, there is a surface 
dead \&yei^, which is much less metallic than the bulk, 
whose thickness C(f7) depends only on bulk properties, 
and diverges for J/buik ^ t^crit in the critical form 

i - (C/crit - C/bulk)"" . (15) 

Therefore £, may be identified with the correlation length 
characteristic of the bulk Mott transition. Numerically, 
we find v = 0.53 ± 0.3 ~ 0.5, a typical mean field critical 
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Figure 4: (Color online) Plot of log(l — -R/-Rbuik) versus z 
for U = 15.97118t (circles), U = 15.9198t (squares), U = 
15.84242t (triangles). In the inset the same data are plotted 
with respect to z (1 — U /UcritY'"^ ■ 



exponent compatible with the simple Gutzwiller approx- 
imation. In Fig. |4] we plot the logarithm of the differ- 
ence between R and i?buik, which clearly shows the ex- 
ponential decay for three values of U. In the inset of the 
same figure we plot the same quantity as function of a 
rescaled coordinate z — *■ z (l — U/UcntY with v = 0.5: all 
data fall on the same curve thus substantiating our state- 
ment on the {/-dependence of the correlation length. Our 
finding of an exponential recovery of the quasiparticle 
weight inside the bulk in place of the expected Friedel- 
like power-law behavior offers a unique opportunity to 
experimentally access the critical properties of the Mott 
transition. Photoemission experiments22 show that the 
surface depletion of metallic electron spectral weight in 
V2O3 propagates inside the interior of the sample for an 
anomalously large depth of many tens of Angstrom be- 
neath the surface, in qualitative agreement with our re- 
sults. Further experiments would be desirable to follow 
the behavior of this length scale upon approaching this 
and other Mott transitions and verify our prediction. 

We end by noting that the calculated Z[z) shows an 
upward curvature near the surface {z = 0), see Fig. [3] 
and also Eq. I|A12P in the appendix. This is unHke earlier 
results obtained by the so-called linearized DMFT— , dis- 
playing instead a linear growth of Z[z) near the surface 
and very close to criticality. Besides a qualitative agree- 
ment with the upward curvature observed in photoemis- 
sion^ which could be coincidental since the real V2O3 is 
much more complicated than our simple one-band Hub- 
bard model, we do not see strong arguments of princi- 
ple supporting either approaches. Both Gutzwiller and 
linearized DMFT are based on rather uncontrolled ap- 
proximations. More reliable techniques, such as straight 
DMFT or Quantum Monte Carlo calculations on large 
size systems, would be needed to clarify this aspect; but 
this is perhaps not important enough. What is more im- 
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Figure 5: Spatial dependence of Z{z) for f/icft = 2t and 
t/right = 15.9712t. The lower panel shows the same data as 
the upper one but closer to the interface. 
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Figure 6: Same as in Fig. [5] for f7icft ~ 15.9198t and C/right 
15.9712£. 



portant is that, just like our approach, also linearized 
DMFT yields, as we checked, to a length controlling the 
depth of the surface perturbed region that diverges at 
the Mott transition. 



B. Geometry (b): Weakly correlated 
metal-strongly correlated metal interface 

The junction between a metal and a Mott insulator or a 
strongly correlated metal was studied recently by Helmes, 
Costi and Roscb^, who used the numerical renormaliza- 
tion group as DMFT impurity solver. With our simpler 
method we can address a broader class of interfaces, in- 
cluding the general case of a correlated metal-correlated 
metal junction, with different values of electron-electron 
interaction in the left (Uioit) and right (J7right) leads. The 



system we consider, see Fig.[ljb), is made of two blocks 
100 layers each, and the junction center is at z = N/2. 
Figs. [Hand [6] show the z dependence of the quasiparticle 
weight for fixed C/right — C4rit and two different values 
of C/icft < C/right- Even if C/(z) is changed stepwise from 
left to right, we find that the closer U\cit is to C/crit, the 
smoother the function Z{z) for z < N/2. On the right 
side of the junction, after a characteristic length bright, 
the quasiparticle weight Z reaches exponentially its bulk 
value. We find for R{z > N/2) a layer dependence well 
represented by the form (for a better fit see Eq. (|A6P with 
the minus sign) 



ht 



left 



i?right)c-(^-^/2)/«'--". (16) 



The dependence of fright on C/right is again given by 
Eq. dm), i.e fright oc (C/^it - C/right)"'' (i^ « 0.5). By 
symmetry, the same holds in the left side too, upon in- 
terchanging the subscripts right and left. 

Our results for weak C/ioft and C/right ^ C/crit can be di- 
rectly compared with those of Helmes et aJ.— , who pro- 
posed that a strongly correlated slab, our right lead with 
C/right — C/crit, in contact with a non interacting metal, 
our left lead, has a quasiparticle weight Z{x) that, close 
to criticality, has a scaling behavior 



x^Z{x) ~ C/ a; 



U - C/crit 



C/, 



crit 



1/2N 



(17) 



where /(O) = 1 and x is the distance from the interface, 
translated in our notation x = z — N/2 and U = C/right- 
The prefactor C ~ 0.008 and the asymptotic behavior 
/(C oo) = 0.15C^ of the scaling function were ex- 
tracted by a DMFT calculation with a 40 layer corre- 
lated slab in contact with a 20 layer almost uncorrelated 
metal^. 

We show in Fig. [7] the quantity Z{x) ex- 
tracted by our Gutzwiller technique and plotted versus 

1 /2 

X |1 — C//C/crit| for different C/'s across the Mott transi- 
tion value. The results are qualitatively similar to those 
of Ref. [33, but differs in two aspects. First of all we find 
that f{() defined in Eq. ifTT)) shows a plateau only when 



2* < a; < 



1 



U 



C/c: 



-1/2 



where an approximate expression for the offset value 
is given in the appendix lA H see Eqs. ljA7P and l|A15p . 
For X <C z*, /(C) ^ so that Z{x) approaches its sur- 
face value at the interface. In our data the crossover 
between the two different regimes is clearly visible, un- 
like in Ref. [s^- More seriously, the coefficient C ~ 0.08 
found by Helmes et aJ.— is almost two orders of magni- 
tude smaller than our, which is numerically around ~ 0.4. 
[The approximate analytical expression discussed in the 
appendix [Al give a slightly larger value of 2/3, see (|Alip 
and I|A17|) ]. In the same appendix we also show that, 
within the Hnearized DMFT approach introduced by Pot- 
thoff and Nolting^S one would extract yet another value 
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Figure 7: (Color online) Plot of Z{x) versus the renor- 
malized coordinate x y^l^^I//C/crit| for U < f/crit (upper 
blue curves: U — 15.7939t triangles, U = 15.84244 crosses, 
U = 15.9198t pluses, U = 15.9712t points, U = 15.9968 tiny 
dots) and U > f/crit (lower blue curves: U = 16.25714 trian- 
gles, U = 16.2035t crosses, U = 16.11484 pluses, U = 16.05114 
points, U — 16.0128 tiny dots). This figure can be compared 
with the inset of Fig. 3 in reference [s^ 



of the coefficient C = 9/11 ^ 0.82, of the same order 
as ours, and again larger than that found by Helmes et 
al^. This disagreement is not just quantitative. Mainly 
because of the smallness of the prefactor, Helmes and 
coworkers^ concluded that the strongly correlated slab 
with U ~ [/crit hence ^buik = .^(x ^ oo) ^ 1 is very 
vifeakly affected by the proximity of the good metal, a 
conclusion later questioned by Zenia et aJ.— , who how- 
ever considered a different geometry. Our results, as well 
as those that could be obtained by linearized DMFT, do 
not allow any such drastic conclusion. Yet, since straight 
DMFT should be more reliable than either linearized 
DMFT or our Gutzwiller approach, it is Hkely that our 
Z{x) is strongly overestimated and that Helmes et aJ.'s 
conclusions are basically correct. It seems worth investi- 
gating further this important question with full DMFT 
on wider slabs. 



C. Geometry (c): Correlated metal-Mott insulator 
(Strongly correlated metal)-correlated metal double 
junction 

In this section we consider geometry (c) of figure [U in 
which a strongly correlated slab of d layers is sandwiched 
between two weakly correlated metal leads, a setup al- 
ready studied by DMFT^>25. In Figs. HI l and [lU] we 
show the layer dependence of the quasiparticle weight for 
different values of the interaction parameters, the Hub- 
bard U in the leads, [/right = f^ieft < t4rit, and in the 



Figure 8: Spatial dependence of Z{z) for [/left = f/right = 24 
and i7ccntcr — 15.97124. The upper panel refers to a central 
region of d = 20 layers, while the lower panel to d = 40 




central slab, Uc, 



> 
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[/crit, and slab thickness d. From 
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Figure 9: Same as in Fig. [H for Uidt = f/right = 15.91984 and 

t/center = 15.97124. 



those results one can draw the following conclusions: 

• For any finite thickness d, the quasiparticle weight 
in the central slab never vanishes, as better revealed 
in Figs. [TT] and [l2l even for [/center > [^crit, fed as it 
is by the evanescent metallic quasiparticle strength 
from the metallic leads. This result agrees perfectly 
with recent DMFT calculations^^. 

• For [/center > [^crit , See Fig. [Tol the minimum value 
Zmin in the central region decreases when d in- 
creases; 

• The behavior of Z{z) across the interface is 
smoother and smoother the closer and closer 

C^right = U\cft are to [/center. 

Looking more in detail at Figs. [H [lO] and at the log- 
scale plots in Fig. [11] and [121 we can identify the char- 
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Figure 10: Same as in Fig. El for (7ieft = Kight = 15.9198* 
and C/ccntor = 16.0288t. 



acteristic differences between a Mott insulating slab and 
a strongly correlated metallic slab, when sandwiched be- 
tween metallic leads. In a strongly correlated metallic 
slab, the central quasiparticle weight ultimately settles 
to the self-standing value it would have in a homoge- 
neous system with U — [/center < C/crit- This value is 
independent of the junction width and of lead correla- 
tions. On the contrary, the quasiparticle weight inside 
the insulating slab is completely borrowed from the leads, 
and strongly depends therefore on their separation and 
correlation. What depends strictly on the central slab in- 
teraction [/center > [/crit is the Quasiparticle decay length 
Ccentor from the lead to the center of the slab, which in- 
creases for increasing slab correlation according to the 
law ([/center — [/crit)"'', with w 0.5, a value that matches 
perfectly that found in section BlI Al 

These considerations suggest that, if we look at the 
problem from a transport point of view, we are con- 
fronted with two completely different mechanisms. In 
a strongly correlated metallic central slab, ^center has the 
role of a screening length, exactly the same role of bright in 
section UlI Bl If instead the central slab is insulating, the 
meaning of ^center becomes completely different, it is now 
a tunneling length. No local quasiparticle peak would 
survive in a homogeneous Mott insulator: the residual 
quasiparticle peak that we find inside the central slab is 
therefore the evanescent lead electron wavefunction that 
tunnels into the slab. 

A special case occurs when [/center ~ [/crit, i-e. right at 
criticality, where neither of the previous two pictures is 
valid. The crossover from the two opposite exponential 
decays describing either screening or tunnehng is char- 
acterized by the absence of any characteristic length, 
which implies a power law variation of the quasiparticle 
strength upon the slab width d 

Z^Ud)-^ + o(^^y (18) 



Figure 11: (Color online) Logarithm of the quasiparticle 
weight Z as St function of layer index z for a 20-sites wide 
(solid line) and 40-sites wide (dashed line) strongly correlated 
metallic slab U = 15.9712t < C/crit sandwiched between two 
weakly correlated metal leads (with U = 15.88438*, 15.79388i, 
15.67674t, 15.53236*.). The entire system is 200-sites wide; 
the interfaces between the leads and the slab are at z = 80 
and z = 120 for the 40-sites wide slab and z = 90 and z = 110 
for the 20-sites wide slab. The figure shows that for increas- 
ing slab width the quasiparticle weight goes to a value that is 
independent of lead correlation. 

We find that the leading l/d^ behavior is, within our 
accuracy, independent of the specific properties of the 
metaUic leads, while the subleading terms do depend on 
them, see Fig. [HI A simple analytical justification of the 
critical behavior is provided in appendix [Al 



IV. CONCLUSIONS 

In this work we have studied how the spatial inho- 
mogeneity of interfaces affects the physics of a strongly 
correlated electron system. To address this problem, 
we extended the conventional Gutzwiller approximation 
technique to account for inhomogeneous Hamiltonian pa- 
rameters. Moreover, to efficiently cope with the larger 
number of variational parameters in comparison with the 
homogeneous case, we derived iterative equations fully 
equivalent to the saddle point equations that identify the 
optimal variational solution, similarly to what is com- 
monly done within unrestricted Hartree-Fock or ab ini- 
tio LDA calculations. These iterative equations can be 
solved without much effort for very large system sizes; an 
advantage with respect to more rigorous approaches, like 
e.g. DMFT calculations, which are numerically feasible 
only for small systems. 

We have applied the method to various interface ge- 
ometries in three dimensions; specifically the interface 
of a strongly correlated metal with the vacuum, the 
interface between two differently correlated metals and 
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Figure 12: (Color online) Same as in Fig. 1111 but the central 
layers have now IJ — 16.1148 > f/crit. In this case the quasi- 
particle weight at the center of the junction is strongly depen- 
dent both on barrier width and on the strength of electron 
correlation in the leads. The central layer remains metallic 
for arbitrary values of C/ > f/crit, but its quasiparticle weight 
decreases exponentially with the slab width. 
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Figure 13: (Color online) Numerical results for and 
V = 15.999t (crosses), 16t (squares), 16.0002t (dashed line), 
16.0004t (diamonds), 16.002t (pluses) for the sandwich geom- 
etry with f7icft = [/right = 2t. The constant value approached 
for U = 16. 00024 ^ Ucrit and large junction width should be 
compared to the one we find in Eq. l(A25|l . 



the junction between two weakly correlated metals sand- 
wiched by a strongly correlated slab. All these geometries 
had been already studied by DMFT21^^>22^^^^^ 
which allowed us to directly compare our results with 
more rigorous ones, thus providing a test on the qual- 
ity of our approximation, which is then applied to much 
larger sizes. 

Our main result is that the effects of an interface de- 
cay exponentially in the interior of a strongly correlated 
system on a very long length-scale proportional to the 
correlation length of the incipient Mott transition, a bulk 
property independent upon the details of the interfacei^ 
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Figure 14: (Color online) Numerical results for Z^ninCp /4: for 
C/icft = [/rigiit ~ 24 [[/center = 16t (squarcs), 16.00024 (dashed 
line), 16.00044 (diamonds)], and for Uidt = [/right = 15.84 
[[/center = 16.00024 (crosses) , 16.00044 (hexagons), 16.00064 
(pluses)]. The stronger lead correlation in the lower curves 
pushes the plateau of the function Zniinrf^/4 towards larger 
values of d. 



In particular, at the surface of a strongly correlated 
metal we find a strong suppression of the metallic prop- 
erties, e.g. of the quasiparticle weight, that persists on a 
large depth controlled by the Mott transition correlation 
length, a "dead layer"— appearing because the surface is 
effectively more correlated than the bulk and consistent 
with photoemission experiments.— Conversely, metallic 
features from a metal lead penetrate inside a Mott in- 
sulator within a depth that, once again, diverges on ap- 
proaching the Mott transition. As a consequence, a con- 
ducting channel always exists inside a Mott insulating 
slab contacted to two metaUic leads, in agreement with 
recent DMFT analyses^^, implying a finite conductance 
at zero bias and temperature that decays fast on increas- 
ing both external parameters on an energy scale expo- 
nentially small in the length of the slab in units of the 
Mott transition correlation length. 

The method that we have developed is very simple and 
flexible, so it can in principle be applied to a variety of re- 
alistic situations of current interest, not only for studying 
interfaces but also for more general inhomogeneities, as 
those arising by impurities or other defects, and can eas- 
ily incorporate additional features like magnetism, which 
we have disregarded throughout this work. 
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Appendix A: ANALYTICAL EXPRESSIONS 
NEAR CRITICALITY 

In this appendix, we show how to derive simple analyt- 
ical expressions for the layer dependence of the quasipar- 
ticle residue near criticality. We assume a three dimen- 
sional slab geometry with constant hopping but inhomo- 



J 



geneous interaction U{z) and with particle-hole symme- 
try. We define as 2e||(z) and 2€±{z — 1/2) the average 
over the uncorrelated Slater determinant |^'o) of the hop- 
ping energy per bond within layer z and between layers 
z and z — 1, respectively. With these definitions, the 
equation ifTTj) can be written as 



2 R{z) (^4: £,1 (z) + e_L(z - 1/2) + e^{z + 1/2)^ + (^e^{z - 1/2) + e^{z + 1/2)^ (^R{z + 1) + R{z - 1) - 2R{z)^ 

U{z) R{z) 



+ e^{z + 1/2) - e^{z ~ 1/2) i?(z + 1) - R{z - 1) + 



r 



^l-R-'iz) 



(Al) 



Near criticality, we expect that the layer dependence 
must appear as a dependence upon the scaling vari- 
able z/^, and, since S, ^ 1, it becomes allowed to re- 
gard z/^ as a continuous variable and expand l|Aip in 
the leading gradients. Because of the interface, both 
e\\{z) and e±[z — 1/2) must acquire a Friedel-like z- 
dependence. However, as shown explicitly in Fig. [HI 
e|| (z) and ej^(z — 1/2) + ej^(z— 1/2) vary appreciably only 
close to the interfaces, while ej_(z — 1/2) — e±{z — 1/2) 
is neghgible. Indeed, as discussed in more detail in the 
Appendix [HI the amplitude of the Friedel's oscillations is 
strongly reduced near criticality, while the period stays 
invariant, so that it is legitimate to neglect the z depen- 
dence of e||(z) and e±{z ± 1/2) and use for them their 
large-z bulk values, en and ej^. 
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Figure 15: Upper panel, plot ekin/t for the sandwich ge- 
ometry (c) with 40 central layers, C/icft = f/right = 2f and 
t/center ~ 15.97124. The value deviates by 2 to 4% from the 
value it would have in a homogeneous system (ekin = t). Mid- 
dle panel, plot of e±sum = e±{z + 1/2) + e±{z — 1/2). Lower 
panel, plot of ej_diff = e±{z + 1/2) — e±{z — 1/2) 



Noting that the average hopping energy per site in the 
homogeneous case is £kin = 4e|| -l-2e_L, the above Eq. (|Aip 
can be written in the continuous limit as 



U 



2R{z)ektn + -r 



R{z) 



4 ^\-R^{2 



2t^ 



d^Rjz) 
dz^ 



0, (A2) 



where we take the bulk value U{z) ~ U, since its vari- 
ation is limited to a single layer. Eq. I|A2|) admits an 
integral of motion, namely 



E 



dR{z) 
dz 



U 



+ Efein R^{z) 

1-^1-R^{z) 

dR{z) 
dz 



E[Riz)]. 



(A3) 



where E [R{z)] is the Gutzwiller variational energy for a 
homogeneous system calculated at fixed R = R{z), i.e. 
not the optimized one. The constant of motion E must 
be chosen to correspond to E[R{zo)] = E[Ro], where zq is 
the layer coordinate at which we expect vanishing deriva- 
tive. In a single interface, we expect that R{z) will reach 
a constant value only asymptotically far from the inter- 
face, i.e. zq ^ CO, where Rq tends to its bulk value 



Ro = ^\-v?, 

and i?[i?o] to the optimized energy in a homogeneous 
system, i.e. 



E = E\Rq\ 



(l-u) 61(1 -u), 



with u = U/Ucrit and Ua-it = ~8^kin, in the Gutzwiller 
approximation. In the case of a correlated slab sand- 
wiched between two metal leads, we expect that R{z) will 
reach a minimum somewhere at midway between the two 
interfaces. If the leads are identical, the minimum occurs 
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right in the middle, so that Rq becomes an unknown pa- 
rameter that has to be fixed by imposing that the actual 
solution R[z,Ro], which depends parametrically on Rq, 
has a vanishing slope dz R[z, Rq] = for z in the middle 
of the slab. 

With the same definitions as above, 



E [Riz)] = + ^ ^ (l - ^Y^Wiz) 

Since in a homogeneous cubic lattice e± = ekin/6 
— t^crit/48, Eq. (|A3P can be rewritten as 



-R^{z) + 2u (l- , (A4) 



where 



Rt, + 2u yjl- = (1 - uY 9(1 - u), (A5) 

in the case of a single interface. The pre-factor 6 of 
the {dR{z) / dz^ comes from the homogeneous relation 
<^kin/^± = 6. As we shall see, the numerical data can be 
better interpreted if Ckin/^i. is considered as a free fitting 
parameter 

The differential equation l|A4| controls the z- 
dependence of R{z > 0), hence of the quasiparticle 
residue Z{z) = R^{z), assuming that the interface affects 
only the boundary condition R{z = 0) = i?surf- There- 
fore, a surface less correlated than the bu lk shou ld be de- 
scribed by l[A4)) with iisurf > i?buik = Vl - u"^ 9(1 - u), 
while by i?surf < -Rbuik the opposite for instance 

the interface with the vacuum of section IIII Al 

We now consider separately the case of a single junc- 
tion and of the double junction, with either metallic or 
insulating bulk. 



1. Single interface with metallic bulk: u < 1 

In the case of a single interface, Eq. IjAsp with u < 1 
has to be used. The differential equation l|A4|) reads 



hence 



namely 



dRiz) 



^ V6 ( y/l^^W{z) - u) , 



= V6i 



dR 



This integral equation can be solved exactly, leading to 
the implicit formula 



V6i 



arc.i„i?,(^) cosa;da; 



arcsin iism-t '"'^^ ^ 

arcsini?(z) — arcsin i?smf 
: tanh"^ 



R{z) i?bull 



tanh 



,i-V(i-^Lk)(i-^^(-)). 

^sm t -Rbulk 



i-\/(i-^buik)(i-^Lt) 



Close to criticality, u ~ 1, one can neglect the arcsines in 
the rhs and find the explicit expression 



Riz) 



i?buik sinhC 



(A6) 



coshC±Vl-^bu 



where the plus sign refers to the case i?surf < -Rbuik, and 
the minus sign to the opposite case, and 



C = ^6(1 - u^) z 



■ tanh 



^surf ^buU 



,i-V(i-<.ik)(i-^Lf) 

= Vg Rbuik {z + z,) . (A7) 

This solution provides a definition of the correlation 
length for u < 1 



1 



v/6(l-ii2) 



0.289 



crit 



C/crit - U 



1/2 



(A8) 



quite close to the DMFT value.^^ We note that, for C > 
1, Eq. I|A6P becomes 



R{z ^ oo) ~ i?buik (^1 T V 1 - ^buik e^'^ 
therefore 

Z{z) = R\z) ^ Zbuik (^1 T 2yr^fl2^e--/« 

■ (A9) 

tends exponentially to its bulk value on a length scale ^, 

from below or above according to i?surf > -Rbuik, respec- 
tively. 

Near criticality, i.e. i?buik = Vl - -C 1, Eq. (|A6l) 
becomes 



R{z) ~ i?buik coth^. 



(AlO) 



so that 



(z + z,)- Z{z) = (z + z,)^ R{zf 
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Figure 16: (Color online) Numerical results for Z{z) in the 
surface geometry, with U = 15.9872t (crosses), 15.9712f (di- 
amonds), 15.9487t (squares), 15.9198t (circles). The solid 
curve is tanh^(C/2), i.e. i?^(C) as defined in Eq. ((Ml) (with 
plus sign) and expanded to first order in -Rbuik ^1. In or- 
der to define the same expansion has been carried out in 
Eq. l|A7p . where we set the quantity Cfein/ex equal to 9.427 
instead of 6, in order to fit the numerical data. 



-(- 

6 V4 



coth^ 



■ /m<i(C): 



(All) 



shows a simple scaling behavior—. The scaHng func- 
tion /u<i(C) that we find has the asymptotic behavior: 
/.<i(0) = 1 and /„<i(C ^ oo) ~ CV4. 



Another case of interest is that of the interface with 
vacuum discussed in section ITlIAl Here i?surf ^ 1 hence 
from Eq. l(A7| it follows that 



■surf 



V6(l -u) 



< 1. 



Away from criticality and for ^ 
since z* <C 1, we find through (jA6[ 
that 



1, which is allowed 
with the plus sign 



so that 



~ a/6 (1 - W) (Z-HZ*): 



Z(z)~6 {\~uf (z + z,)^ 



(A12) 



showing that the quasiparticle residue approaches its sur- 
face value with a finite curvature. 

In Fig. [inland Fig.[l7|we show that rescaled numerical 
data for an interface between a 200-layer-wide correlated 
metal slab and the vacuum and for a junction between a 
weakly correlated metal and a strongly correlated metal. 



It is easy to fit the numerical data with the function 
i?^(z) displayed in Eq. I|A7P by tuning just one param- 
eter, which, as discussed above, is the value of tkml^i. 
(equal to 6 in the homogeneous problem). The fact that 
the ideal theoretical result, relying on homogeneous val- 
ues for hopping and kinetic energy, fits the numerical 
data with just a single tunable parameter, is a pleasant 
feature. 




Figure 17: (Color online) Numerical results for Z(^z) in the 
single junction geometry with metallic bulk, the position of 
the junction is chosen as the origin for the spatial coordinate, 
the metal on the left side is very weakly correlated (C/ — 2t); 
the values for U on the right side are the same of Fig. 1161 The 
solid curve is now the function 1/ tanh^((^/2), i.e. the second 



power of Eq. IIA6P (with minus sign) expanded to first order 
in -Rbuik. As in Fig. 1161 the definition of has been obtained 
from Eq. I|A7P by expanding to first order in 7?buik. The value 
of Efcin/ex that fits the data is now 8.254. 



2. Single interface with insulating bulk: u > 1 

In this case the equation l|A4[) using l|A5[) with w > 1 

reads 



1 ( dR{z) 
leading to 



(A13) 



/ 



dR 



2u~R^ - 2u\/l - i?2 



= — \/6 J dz, 



where we have assumed that on the surface i?surf is finite 
and decay in the bulk, so that the derivative is negative. 
The above integral equation can be solved too, with an 
implicit solution 
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-•\/6 (w — I) z — 2Vu — 1 arcsin ( '^"^ ^""^ ) — 2 Vm - 



tanh 



— 1 cos 
\J u — cos^ 

I 



tanh 



1 arcsin 



cos J/surf 



1 cos J/surf 



COS^ ysurf 



where = sin2j/(z), i?surf = sm2ysurf- As before the 
arcsines can be neglected near criticahty to obtain the 
explicit solution 



R\z) = 1 - 1 - 



2(w- 1) 
u cosh^ C ~ 1 / 



(A14) 



with 



_1 / - 1 COSysurt 



C = a/B (u — 1) z + tanh 

COS^ J/surf / 

= v/6(w- 1) (z + z*). (A15) 

In the case of an insulating bulk, the correlation length 
defined through (|A16P is therefore 



^6 (m- 1) 



0.408 



C/c: 



1/2 



(A16) 



with a different numerical prefactor, actually a -\/2 
greater, with respect to the metahic bulk IjAsp . 
Near criticality, w > 1, 



= Z{z) 



4(u- 1) 
sinh^C ' 



so that, as before. 



j/.mK), (A17) 



has a scaling behavior with /u>i(0) = 1 and 
/n>i(C^oo)~4C'e-2^. 

3. Double junction 

We assume for simplicity a slab of length 2L in contact 
with two leads. In this case we need to use Eq. (jA4p 
with i?o a parameter that has to be fixed by imposing 
that the solution R{z) becomes i?o at some zq within the 
slab. If we assume that both leads are less correlated 
than the slab, then R(z) always decreases moving away 
from any of the two interfaces, and we can determine i?o 
by imposing either of the two fohowing conditions: 



dR 



Y ^0 + 2u-y/l - Rl - i?2 - 2uVl - i?2 
dR 

'rI + 2u^l -Rl-R^^ 2uVl - R^ 

I 



= -VGzo, (A18) 
= V6(2L-zo),, (A19) 



where -R^^.f and R^^^f are the values of R{z) at the left (|A19P minus (|A18P we find 
and right surfaces, respectively. Taking the difference 



V62L : 




dR 



Rl + 2uxll 



(A20) 



/?2 



i?2 - 2u\/l - i?2 
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which has to be solved to find Rq as function of the other 
parameters. Once Rq is found, one can determine zg. In 
order to simpHfy the calculations, we will assume two 



identical leads, i.e. R^^^f = ^^,i.f = -Rsurf , so that zq = L 
and l|A20p becomes 



V6L 



dR 



R', 



2uJl - Rl- R^ - 2u\/l - i?2 



v/(a-c)(6-rf) 



(c - 6) n 



b-d 



bFi^,k) 



iA211 



with parameters a>b>c>u>d. The last ex- 
pression can be de rived easily after the cha nge of vari- 
able R = Vl — x'^, and seemingly Rq = — x\ and 
i?suit = v^l — a:^surf n(0, n, fc) and F{(l),k) are elliptic 
integrals of third and first kind, respectively 



Tl{4>, n, k) 



dx 



- fc^ sin^ X 
dx 



(1 



n sm X 



sin^ X 



and 



'(6-d)(c-u) 
{c-d){b-uy 



k = 



'{a-b){c~d) 



{a-c){b-d) 
The various parameters are, when 2u — xq > 1, 

a = 2u — xq, 

b = 1, 

c = xo, 

d = -1, 

W -^surf ; 

so that 



arcsm < 



2 (xq ~ Xsurf) 
(xo + 1) (1 - Xsurf ) ' 



'(2M-a;o-l)(a;o + l) 



y 4 (it - Xo) 

On the contrary, if 2m — xo < 1 , then 



a ~ 1, 

b = 2u — xo, 

c = Xo, 



d = -1, 



hence 



(2u -Xq + 1){xo - Xsurl) 
{xq + l){2u - Xo - Xsurf) '' 



We rewrite 
(c - 5) n 

dx 



'(1-2u + .to)(xo + 1) 
(1 - xo){2u~xo + 1)' 



+6F((/.,fc) (A22) 

d{b - c) + b{c - d) cos^ x\ 1 
{b-c) + {c-d)cof?x J Vl - fc2 sin^x' 



and note that Sit x = 4> 

d{b - c) + 5(c - d) cos^ 
(6 - c) + (c - d) cos^ (j) 



In addition 6 — c in both cases is very small. Indeed, 
for 2u — Xo > 1, whic h corresponds to an insulating slab 
where Rq = y/l — x^ for large L, b ~ c — 1 ~ xo -^i 
1. In the opposite case of a weakly correlated slab, still 
6 — c = 2u — Xo — Xo ^ 1 since xq ^ u for large L. 
Therefore 

d{b — c) + b{c — d) cos^ x 
(6 — c) + (c — d) cos^ X 

is practically constant and equal to b everywhere but 
close to the extreme of integration, where it fastly de- 
cays to Xsurf- Therefore to leading order we can write 



(c - 6) n 



b-d' 



,k] +bF{(f>,k) ~ bF{(j),k), 



hence the equation to be solved becomes 
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V6L 



2b 



^{a~c){b-d) 



2b 



V{a-c){b-d) 



K{k) — F arcsin 



cos ( 



1 — /c^ sin^ I 



(A23) 



where K{k) = F{Tr/2,k) and the last expression being 
more convenient since (j) ~ 7r/2. 

In order to find xq as function of the other parameters, 
we have to consider separately three different cases. 



a. Insulating off-critical behavior: u 2> 1 

In this case 2u — xq > 1 . We note that fc as a function 
of u at fixed xo — 1 is equal to 

;,2 _ + 1 1 



for u = 1, and very rapidly increases to its asymptotic 
M ^ 1 value 



Therefore (|A23P is, at leading order. 



K 



I+Xq 



v^:^ \\ 2 

Therefore, in this limit. 



1 , 32 
m ■ 



2Vu - 1 1 - a;o ' 



Zo = Ro^ 64 c~V24(«-i) (^24) 
vanishes exponentially in the length of the slab. 



In this case 



b. Critical behavior: it = 1 



4 "2' 

hence at leading order Eq. (|A23P reads 



V6L = 



1 



K 



1 



1 



Vl — 2^0 W2/ 4\/7r Vl — Xq 
from which it follows that 
1 



Zq — Rq ^ 



487r 



1 4 



1 

l2 



1.146 



(A25) 



Once again we find a critical behavior L^Zq ~ const., 
with a sizable constant 1.146. 



c. Metallic off-critical behavior: u <C 1 

This is the case in which 2u — xq < 1 and xq ~ m, so 
that 



1-k'^ 



4(u - Xq) 



4(u - xo) 



(2u- xi + 1)(1 -xo) 1 



.2 ■ 



Therefore Eq. l(A23l) is 



In- 



u ^^4(l-u2) 



a/1 -u2 1-A:^ Vl -u^ 



whose solution is 

M - Xo = 4(1 - c-v^v^T=^ V". 

Therefore, since Zbuik = 1 — u^, it follows that 

Zo ~ Zbuik (l + 8w e-^^/^^/") . (A26) 

4. Comparison with DMFT 

Near the Mott transition, U ~ ?7crit, Potthoff and Nolt- 
ing in Ref. [H have introduced a set of linearized DMFT 
recursive equations for the layer dependent quasiparticle 
residue. Taking, as before, the continuous limit of their 
Eq. (37), with q = A p = 1 and C/crit — 6t-\/6, one finds 
the following differential equation 



1 d^Zjz) 
6 dz^ 



+ 2Z{z) {\-u)-cZ{zf ^Q. (A27) 



The numerical constant is estimated to be c = 11/9^. 
The limiting behavior for u ^ 1 is the solution of 



1 d'^z{z) 

6 dz^ 



cZ{z)\ 



namely 



9 , ^ 1 9 
.^Z(.)^- = -.0.82. 



(A28) 



Let's consider instead our Eq. \K2} that, divided by 

2tkin = — t^crit/4. Can bc written as 

1 d'^R{z) ^, ^ R{z) 
= - „ T + R{z) - u ■ ' 



6 dz' 
_ 1 d'^Ri^z) 
~ 6 dz'^ 



+ (l-u) R{z)--R{zf.{k2<=>) 
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Figure 18: (Color online) Mean value of the hopping matrix 
element on the uncorrelated wavefunction versus the distance 
from the leftmost surface layer in geometry (a) with C/surfaco ~ 
20t and (7buik = 14.6642t (triangles, panel 1) and (7buik = 
15.9712t (squares, panel 3). The circles in panel 2 show the 
hopping for the same simulation that was performed for panel 
1, but as a function of distance from the right surface of the 
sample, where U = t/buik = 14.6642t. The results of fit are 
showed by the solid lines. From above, the first and second 
curves are a plot of Eq. ([Ml with A = 0.1673, w = -0.0046 
and A = 0.1673, w = —0.0074 respectively. 



At criticality, u ^ 1, the solution 



z^Z{z) - ~ 0.66, 
o 



(A30) 



is just the Hmiting value of Eqs. (|A17P and (jAlip for 
( = 0. The numerical coefficient 2/3 that we find is 
sHghtly smaller than the linearized DMFT one, 9/11, 
but both are much bigger than the value extracted by 
straight DMFT calculations in Ref. [H, namely 0.008. 
Supposedly, straight DMFT is a better approximation 
than Hnearized DMFT, which in turns should be better 
than our Gutzwiller technique, therefore it is likely that 
our results overestimate the quasiparticle residue Z. 



Appendix B: FRIEDEL'S OSCILLATIONS 

In the previous sections we have derived a simple model 
to extract the behavior of Z{z) assuming uniform val- 
ues for the hopping matrix elements on the uncorrelated 
Slater determinant. Of course the hopping is not uni- 
form, its variation being described in most cases by some 
Friedel oscillations around the bulk value (thin solid lines 
in Figs. [2l fT0|) . The Friedel's oscillations arise as a conse- 
quence of broken translational symmetry in a Fermi gas, 
i.e. around a single impurity or near an interface. An 



impurity embedded in an electron gas of dimensionality 
D induces oscillations that decay as a power law 
and whose wavevector is twice the Fermi wavevector— . 
The Friedel's oscillations in a D = 3 electron gas with 
an interface can be obtained as a superposition of Friedel 
oscillations for a layer of impurities, and one can readily 
find that, moving perpendicularly to the interface over a 
length X, they behave at leading order as 



cos 2kpx 



(Bl) 



results which is strictly valid for a spherical Fermi surface, 
although the decay exponent is independent of the shape 
of the Fermi surface. 

If we include electron-electron interaction via the Hub- 
bard U and treat it by the Gutzwiller approximation, we 
expect that the Friedel's oscillation will be affected also 
by the layer-dependence of the quasiparticle weight Z{z). 
Our results show that the faster the change of Z{z), the 
larger the oscillations. This means that a system with 
geometry (a) and C/buik ~ C^crit displays much smoother 
oscillations that a system with J/buik ^ C4rit, since the 
spatial dependence of Z{z) is sharper when the bulk in- 
teraction parameter is far from criticality. 

In light of the spatial dependence of the oscillations 
predicted by Eq. l|Bip . we fitted our data for the hop- 
ping e±{x + 1/2) perpendicular to the interface and in 
geometries (a) and (b) (see Fig. [l]) with the function 



A 



cos TTX 



(B2) 



where x is the distance from either the surface layer (ge- 
ometry (a)) or the layer across which U{z) changes step- 
wise (geometry (b)). The function (jB2p fits the data 
showed in Figs. [18] for a weakly correlated system with 
strongly correlated surface. If the bulk value of U is in- 
creased towards J7crit , the correlation length ^ becomes so 
big that it is hard to identify unambiguously any Friedel's 
oscillation, as shown in Fig. [18] panel (3). The function 
(|B2p fits also the data for the hopping on the weakly cor- 
related side of the junction in geometry (b), see Fig. [T9|) . 
On the strongly correlated (right) side of the same junc- 
tion again the correlation length ^ is too large and we 
were not able to make any fit. 

In conclusion, the inhomogeneity of the interaction pa- 
rameter U affects the spatial dependence not only of the 
quasiparticle weight, but also of the hopping matrix ele- 
ment on the uncorrelated Slater determinant. The latter 
displays Friedel's oscillations that rise from the breaking 
of discrete translational symmetry. In any case, when 
the system is in the close vicinity of the Mott transition, 
the effects of these oscillations are smoothed out as a re- 
sult of the diverging characteristic length ^ of the local 
quasiparticle weight. 
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Figure 19: (Color online) Plot of the hopping matrix ele- 
ment for a system with geometry (b), C/icft = 2t and (/right = 
15.97118t. Upper panel: x is the distance from the junction 
on the weakly correlated metallic (left) side; lower panel: the 
same on the strongly correlated metallic (right) side. The 
Friedel oscillations on the weakly correlated side are fitted by 
Eq. (dU with A = 0.16715, w = 0.0050. On the strongly cor- 
related side the fit was not possible for the reasons explained 
in the text. 
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